/*
 *   This program is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 *    PaceMatrix.java
 *    Copyright (C) 2002 University of Waikato, Hamilton, New Zealand
 *
 */

package weka.classifiers.functions.pace;

import java.text.DecimalFormat;
import java.util.Random;

import weka.core.RevisionUtils;
import weka.core.matrix.DoubleVector;
import weka.core.matrix.FlexibleDecimalFormat;
import weka.core.matrix.IntVector;
import weka.core.matrix.Maths;
import weka.core.matrix.Matrix;

/**
 * Class for matrix manipulation used for pace regression.
 * <p>
 * 
 * REFERENCES
 * <p>
 * 
 * Wang, Y. (2000). "A new approach to fitting linear models in high dimensional
 * spaces." PhD Thesis. Department of Computer Science, University of Waikato,
 * New Zealand.
 * <p>
 * 
 * Wang, Y. and Witten, I. H. (2002). "Modeling for optimal probability
 * prediction." Proceedings of ICML'2002. Sydney.
 * <p>
 * 
 * @author Yong Wang (yongwang@cs.waikato.ac.nz)
 * @version $Revision$
 */
public class PaceMatrix extends Matrix {

  /** for serialization */
  static final long serialVersionUID = 2699925616857843973L;

  /*
   * ------------------------ Constructors ------------------------
   */

  /**
   * Construct an m-by-n PACE matrix of zeros.
   * 
   * @param m Number of rows.
   * @param n Number of colums.
   */
  public PaceMatrix(int m, int n) {
    super(m, n);
  }

  /**
   * Construct an m-by-n constant PACE matrix.
   * 
   * @param m Number of rows.
   * @param n Number of colums.
   * @param s Fill the matrix with this scalar value.
   */
  public PaceMatrix(int m, int n, double s) {
    super(m, n, s);
  }

  /**
   * Construct a PACE matrix from a 2-D array.
   * 
   * @param A Two-dimensional array of doubles.
   * @throws IllegalArgumentException All rows must have the same length
   */
  public PaceMatrix(double[][] A) {
    super(A);
  }

  /**
   * Construct a PACE matrix quickly without checking arguments.
   * 
   * @param A Two-dimensional array of doubles.
   * @param m Number of rows.
   * @param n Number of colums.
   */
  public PaceMatrix(double[][] A, int m, int n) {
    super(A, m, n);
  }

  /**
   * Construct a PaceMatrix from a one-dimensional packed array
   * 
   * @param vals One-dimensional array of doubles, packed by columns (ala
   *          Fortran).
   * @param m Number of rows.
   * @throws IllegalArgumentException Array length must be a multiple of m.
   */
  public PaceMatrix(double vals[], int m) {
    super(vals, m);
  }

  /**
   * Construct a PaceMatrix with a single column from a DoubleVector
   * 
   * @param v DoubleVector
   */
  public PaceMatrix(DoubleVector v) {
    this(v.size(), 1);
    setMatrix(0, v.size() - 1, 0, v);
  }

  /**
   * Construct a PaceMatrix from a Matrix
   * 
   * @param X Matrix
   */
  public PaceMatrix(Matrix X) {
    super(X.getRowDimension(), X.getColumnDimension());
    A = X.getArray();
  }

  /*
   * ------------------------ Public Methods ------------------------
   */

  /**
   * Set the row dimenion of the matrix
   * 
   * @param rowDimension the row dimension
   */
  public void setRowDimension(int rowDimension) {
    m = rowDimension;
  }

  /**
   * Set the column dimenion of the matrix
   * 
   * @param columnDimension the column dimension
   */
  public void setColumnDimension(int columnDimension) {
    n = columnDimension;
  }

  /**
   * Clone the PaceMatrix object.
   * 
   * @return the clone
   */
  @Override
  public Object clone() {
    PaceMatrix X = new PaceMatrix(m, n);
    double[][] C = X.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        C[i][j] = A[i][j];
      }
    }
    return X;
  }

  /**
   * Add a value to an element and reset the element
   * 
   * @param i the row number of the element
   * @param j the column number of the element
   * @param s the double value to be added with
   */
  public void setPlus(int i, int j, double s) {
    A[i][j] += s;
  }

  /**
   * Multiply a value with an element and reset the element
   * 
   * @param i the row number of the element
   * @param j the column number of the element
   * @param s the double value to be multiplied with
   */
  public void setTimes(int i, int j, double s) {
    A[i][j] *= s;
  }

  /**
   * Set the submatrix A[i0:i1][j0:j1] with a same value
   * 
   * @param i0 the index of the first element of the column
   * @param i1 the index of the last element of the column
   * @param j0 the index of the first column
   * @param j1 the index of the last column
   * @param s the value to be set to
   */
  public void setMatrix(int i0, int i1, int j0, int j1, double s) {
    try {
      for (int i = i0; i <= i1; i++) {
        for (int j = j0; j <= j1; j++) {
          A[i][j] = s;
        }
      }
    } catch (ArrayIndexOutOfBoundsException e) {
      throw new ArrayIndexOutOfBoundsException("Index out of bounds");
    }
  }

  /**
   * Set the submatrix A[i0:i1][j] with the values stored in a DoubleVector
   * 
   * @param i0 the index of the first element of the column
   * @param i1 the index of the last element of the column
   * @param j the index of the column
   * @param v the vector that stores the values
   */
  public void setMatrix(int i0, int i1, int j, DoubleVector v) {
    for (int i = i0; i <= i1; i++) {
      A[i][j] = v.get(i - i0);
    }
  }

  /**
   * Set the whole matrix from a 1-D array
   * 
   * @param v 1-D array of doubles
   * @param columnFirst Whether to fill the column first or the row.
   * @throws ArrayIndexOutOfBoundsException Submatrix indices
   */
  public void setMatrix(double[] v, boolean columnFirst) {
    try {
      if (v.length != m * n) {
        throw new IllegalArgumentException("sizes not match.");
      }
      int i, j, count = 0;
      if (columnFirst) {
        for (i = 0; i < m; i++) {
          for (j = 0; j < n; j++) {
            A[i][j] = v[count];
            count++;
          }
        }
      } else {
        for (j = 0; j < n; j++) {
          for (i = 0; i < m; i++) {
            A[i][j] = v[count];
            count++;
          }
        }
      }

    } catch (ArrayIndexOutOfBoundsException e) {
      throw new ArrayIndexOutOfBoundsException("Submatrix indices");
    }
  }

  /**
   * Returns the maximum absolute value of all elements
   * 
   * @return the maximum value
   */
  public double maxAbs() {
    double ma = Math.abs(A[0][0]);
    for (int j = 0; j < n; j++) {
      for (int i = 0; i < m; i++) {
        ma = Math.max(ma, Math.abs(A[i][j]));
      }
    }
    return ma;
  }

  /**
   * Returns the maximum absolute value of some elements of a column, that is,
   * the elements of A[i0:i1][j].
   * 
   * @param i0 the index of the first element of the column
   * @param i1 the index of the last element of the column
   * @param j the index of the column
   * @return the maximum value
   */
  public double maxAbs(int i0, int i1, int j) {
    double m = Math.abs(A[i0][j]);
    for (int i = i0 + 1; i <= i1; i++) {
      m = Math.max(m, Math.abs(A[i][j]));
    }
    return m;
  }

  /**
   * Returns the minimum absolute value of some elements of a column, that is,
   * the elements of A[i0:i1][j].
   * 
   * @param i0 the index of the first element of the column
   * @param i1 the index of the last element of the column
   * @param column the index of the column
   * @return the minimum value
   */
  public double minAbs(int i0, int i1, int column) {
    double m = Math.abs(A[i0][column]);
    for (int i = i0 + 1; i <= i1; i++) {
      m = Math.min(m, Math.abs(A[i][column]));
    }
    return m;
  }

  /**
   * Check if the matrix is empty
   * 
   * @return true if the matrix is empty
   */
  public boolean isEmpty() {
    if (m == 0 || n == 0) {
      return true;
    }
    if (A == null) {
      return true;
    }
    return false;
  }

  /**
   * Return a DoubleVector that stores a column of the matrix
   * 
   * @param j the index of the column
   * @return the column
   */
  public DoubleVector getColumn(int j) {
    DoubleVector v = new DoubleVector(m);
    double[] a = v.getArray();
    for (int i = 0; i < m; i++) {
      a[i] = A[i][j];
    }
    return v;
  }

  /**
   * Return a DoubleVector that stores some elements of a column of the matrix
   * 
   * @param i0 the index of the first element of the column
   * @param i1 the index of the last element of the column
   * @param j the index of the column
   * @return the DoubleVector
   */
  public DoubleVector getColumn(int i0, int i1, int j) {
    DoubleVector v = new DoubleVector(i1 - i0 + 1);
    double[] a = v.getArray();
    int count = 0;
    for (int i = i0; i <= i1; i++) {
      a[count] = A[i][j];
      count++;
    }
    return v;
  }

  /**
   * Multiplication between a row (or part of a row) of the first matrix and a
   * column (or part or a column) of the second matrix.
   * 
   * @param i the index of the row in the first matrix
   * @param j0 the index of the first column in the first matrix
   * @param j1 the index of the last column in the first matrix
   * @param B the second matrix
   * @param l the index of the column in the second matrix
   * @return the result of the multiplication
   */
  public double times(int i, int j0, int j1, PaceMatrix B, int l) {
    double s = 0.0;
    for (int j = j0; j <= j1; j++) {
      s += A[i][j] * B.A[j][l];
    }
    return s;
  }

  /**
   * Decimal format for converting a matrix into a string
   * 
   * @return the default decimal format
   */
  protected DecimalFormat[] format() {
    return format(0, m - 1, 0, n - 1, 7, false);
  }

  /**
   * Decimal format for converting a matrix into a string
   * 
   * @param digits the number of digits
   * @return the decimal format
   */
  protected DecimalFormat[] format(int digits) {
    return format(0, m - 1, 0, n - 1, digits, false);
  }

  /**
   * Decimal format for converting a matrix into a string
   * 
   * @param digits the number of digits
   * @param trailing
   * @return the decimal format
   */
  protected DecimalFormat[] format(int digits, boolean trailing) {
    return format(0, m - 1, 0, n - 1, digits, trailing);
  }

  /**
   * Decimal format for converting a matrix into a string
   * 
   * @param i0
   * @param i1
   * @param j
   * @param digits the number of digits
   * @param trailing
   * @return the decimal format
   */
  protected DecimalFormat format(int i0, int i1, int j, int digits,
    boolean trailing) {
    FlexibleDecimalFormat df = new FlexibleDecimalFormat(digits, trailing);
    df.grouping(true);
    for (int i = i0; i <= i1; i++) {
      df.update(A[i][j]);
    }
    return df;
  }

  /**
   * Decimal format for converting a matrix into a string
   * 
   * @param i0
   * @param i1
   * @param j0
   * @param j1
   * @param trailing
   * @param digits the number of digits
   * @return the decimal format
   */
  protected DecimalFormat[] format(int i0, int i1, int j0, int j1, int digits,
    boolean trailing) {
    DecimalFormat[] f = new DecimalFormat[j1 - j0 + 1];
    for (int j = j0; j <= j1; j++) {
      f[j] = format(i0, i1, j, digits, trailing);
    }
    return f;
  }

  /**
   * Converts matrix to string
   * 
   * @return the matrix as string
   */
  @Override
  public String toString() {
    return toString(5, false);
  }

  /**
   * Converts matrix to string
   * 
   * @param digits number of digits after decimal point
   * @param trailing true if trailing zeros are padded
   * @return the matrix as string
   */
  public String toString(int digits, boolean trailing) {

    if (isEmpty()) {
      return "null matrix";
    }

    StringBuffer text = new StringBuffer();
    DecimalFormat[] nf = format(digits, trailing);
    int numCols = 0;
    int count = 0;
    int width = 80;
    int lenNumber;

    int[] nCols = new int[n];
    int nk = 0;
    for (int j = 0; j < n; j++) {
      lenNumber = nf[j].format(A[0][j]).length();
      if (count + 1 + lenNumber > width - 1) {
        nCols[nk++] = numCols;
        count = 0;
        numCols = 0;
      }
      count += 1 + lenNumber;
      ++numCols;
    }
    nCols[nk] = numCols;

    nk = 0;
    for (int k = 0; k < n;) {
      for (int i = 0; i < m; i++) {
        for (int j = k; j < k + nCols[nk]; j++) {
          text.append(" " + nf[j].format(A[i][j]));
        }
        text.append("\n");
      }
      k += nCols[nk];
      ++nk;
      text.append("\n");
    }

    return text.toString();
  }

  /**
   * Squared sum of a column or row in a matrix
   * 
   * @param j the index of the column or row
   * @param i0 the index of the first element
   * @param i1 the index of the last element
   * @param col if true, sum over a column; otherwise, over a row
   * @return the squared sum
   */
  public double sum2(int j, int i0, int i1, boolean col) {
    double s2 = 0;
    if (col) { // column
      for (int i = i0; i <= i1; i++) {
        s2 += A[i][j] * A[i][j];
      }
    } else {
      for (int i = i0; i <= i1; i++) {
        s2 += A[j][i] * A[j][i];
      }
    }
    return s2;
  }

  /**
   * Squared sum of columns or rows of a matrix
   * 
   * @param col if true, sum over columns; otherwise, over rows
   * @return the squared sum
   */
  public double[] sum2(boolean col) {
    int l = col ? n : m;
    int p = col ? m : n;
    double[] s2 = new double[l];
    for (int i = 0; i < l; i++) {
      s2[i] = sum2(i, 0, p - 1, col);
    }
    return s2;
  }

  /**
   * Constructs single Householder transformation for a column
   * 
   * @param j the index of the column
   * @param k the index of the row
   * @return d and q
   */
  public double[] h1(int j, int k) {
    double dq[] = new double[2];
    double s2 = sum2(j, k, m - 1, true);
    dq[0] = A[k][j] >= 0 ? -Math.sqrt(s2) : Math.sqrt(s2);
    A[k][j] -= dq[0];
    dq[1] = A[k][j] * dq[0];
    return dq;
  }

  /**
   * Performs single Householder transformation on one column of a matrix
   * 
   * @param j the index of the column
   * @param k the index of the row
   * @param q q = - u'u/2; must be negative
   * @param b the matrix to be transformed
   * @param l the column of the matrix b
   */
  public void h2(int j, int k, double q, PaceMatrix b, int l) {
    double s = 0, alpha;
    for (int i = k; i < m; i++) {
      s += A[i][j] * b.A[i][l];
    }
    alpha = s / q;
    for (int i = k; i < m; i++) {
      b.A[i][l] += alpha * A[i][j];
    }
  }

  /**
   * Constructs the Givens rotation
   * 
   * @param a
   * @param b
   * @return a double array that stores the cosine and sine values
   */
  public double[] g1(double a, double b) {
    double cs[] = new double[2];
    double r = Maths.hypot(a, b);
    if (r == 0.0) {
      cs[0] = 1;
      cs[1] = 0;
    } else {
      cs[0] = a / r;
      cs[1] = b / r;
    }
    return cs;
  }

  /**
   * Performs the Givens rotation
   * 
   * @param cs a array storing the cosine and sine values
   * @param i0 the index of the row of the first element
   * @param i1 the index of the row of the second element
   * @param j the index of the column
   */
  public void g2(double cs[], int i0, int i1, int j) {
    double w = cs[0] * A[i0][j] + cs[1] * A[i1][j];
    A[i1][j] = -cs[1] * A[i0][j] + cs[0] * A[i1][j];
    A[i0][j] = w;
  }

  /**
   * Forward ordering of columns in terms of response explanation. On input,
   * matrices A and b are already QR-transformed. The indices of transformed
   * columns are stored in the pivoting vector.
   * 
   * @param b the PaceMatrix b
   * @param pvt the pivoting vector
   * @param k0 the first k0 columns (in pvt) of A are not to be changed
   **/
  public void forward(PaceMatrix b, IntVector pvt, int k0) {
    for (int j = k0; j < Math.min(pvt.size(), m); j++) {
      steplsqr(b, pvt, j, mostExplainingColumn(b, pvt, j), true);
    }
  }

  /**
   * Returns the index of the column that has the largest (squared) response,
   * when each of columns pvt[ks:] is moved to become the ks-th column. On
   * input, A and b are both QR-transformed.
   * 
   * @param b response
   * @param pvt pivoting index of A
   * @param ks columns pvt[ks:] of A are to be tested
   * @return the index of the column
   */
  public int mostExplainingColumn(PaceMatrix b, IntVector pvt, int ks) {
    double val;
    // int [] p = pvt.getArray(); NOT USED
    double ma = columnResponseExplanation(b, pvt, ks, ks);
    int jma = ks;
    for (int i = ks + 1; i < pvt.size(); i++) {
      val = columnResponseExplanation(b, pvt, i, ks);
      if (val > ma) {
        ma = val;
        jma = i;
      }
    }
    return jma;
  }

  /**
   * Backward ordering of columns in terms of response explanation. On input,
   * matrices A and b are already QR-transformed. The indices of transformed
   * columns are stored in the pivoting vector.
   * 
   * A and b must have the same number of rows, being the (pseudo-)rank.
   * 
   * @param b PaceMatrix b
   * @param pvt pivoting vector
   * @param ks number of QR-transformed columns; psuedo-rank of A
   * @param k0 first k0 columns in pvt[] are not to be ordered.
   */
  public void backward(PaceMatrix b, IntVector pvt, int ks, int k0) {
    for (int j = ks; j > k0; j--) {
      steplsqr(b, pvt, j, leastExplainingColumn(b, pvt, j, k0), false);
    }
  }

  /**
   * Returns the index of the column that has the smallest (squared) response,
   * when the column is moved to become the (ks-1)-th column. On input, A and b
   * are both QR-transformed.
   * 
   * @param b response
   * @param pvt pivoting index of A
   * @param ks psudo-rank of A
   * @param k0 A[][pvt[0:(k0-1)]] are excluded from the testing.
   * @return the index of the column
   */
  public int leastExplainingColumn(PaceMatrix b, IntVector pvt, int ks, int k0) {
    double val;
    // int [] p = pvt.getArray(); NOT USED
    double mi = columnResponseExplanation(b, pvt, ks - 1, ks);
    int jmi = ks - 1;
    for (int i = k0; i < ks - 1; i++) {
      val = columnResponseExplanation(b, pvt, i, ks);
      if (val <= mi) {
        mi = val;
        jmi = i;
      }
    }
    return jmi;
  }

  /**
   * Returns the squared ks-th response value if the j-th column becomes the
   * ks-th after orthogonal transformation. A[][pvt[ks:j]] (or A[][pvt[j:ks]],
   * if ks > j) and b[] are already QR-transformed on input and will remain
   * unchanged on output.
   * 
   * More generally, it returns the inner product of the corresponding row
   * vector of the response PaceMatrix. (To be implemented.)
   * 
   * @param b PaceMatrix b
   * @param pvt pivoting vector
   * @param j the column A[pvt[j]][] is to be moved
   * @param ks the target column A[pvt[ks]][]
   * @return the squared response value
   */
  public double columnResponseExplanation(PaceMatrix b, IntVector pvt, int j,
    int ks) {
    /*
     * Implementation:
     * 
     * If j == ks - 1, returns the squared ks-th response directly.
     * 
     * If j > ks -1, returns the ks-th response after Householder-transforming
     * the j-th column and the response.
     * 
     * If j < ks - 1, returns the ks-th response after a sequence of Givens
     * rotations starting from the j-th row.
     */

    int k, l;
    double[] xxx = new double[n];
    int[] p = pvt.getArray();
    double val;

    if (j == ks - 1) {
      val = b.A[j][0];
    } else if (j > ks - 1) {
      int jm = Math.min(n - 1, j);
      DoubleVector u = getColumn(ks, jm, p[j]);
      DoubleVector v = b.getColumn(ks, jm, 0);
      val = v.innerProduct(u) / u.norm2();
    } else { // ks > j
      for (k = j + 1; k < ks; k++) {
        xxx[k] = A[j][p[k]];
      }
      val = b.A[j][0];
      double[] cs;
      for (k = j + 1; k < ks; k++) {
        cs = g1(xxx[k], A[k][p[k]]);
        for (l = k + 1; l < ks; l++) {
          xxx[l] = -cs[1] * xxx[l] + cs[0] * A[k][p[l]];
        }
        val = -cs[1] * val + cs[0] * b.A[k][0];
      }
    }
    return val * val; // or inner product in later implementation???
  }

  /**
   * QR transformation for a least squares problem<br/>
   * A x = b<br/>
   * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of
   * A.
   * 
   * @param b PaceMatrix b
   * @param pvt pivoting vector
   * @param k0 the first k0 columns of A (indexed by pvt) are pre-chosen. (But
   *          subject to rank examination.)
   * 
   *          For example, the constant term may be reserved, in which case k0 =
   *          1.
   **/
  public void lsqr(PaceMatrix b, IntVector pvt, int k0) {
    final double TINY = 1e-15;
    int[] p = pvt.getArray();
    int ks = 0; // psuedo-rank
    for (int j = 0; j < k0; j++) {
      if (sum2(p[j], ks, m - 1, true) > TINY) { // large diagonal element
        steplsqr(b, pvt, ks, j, true);
        ks++;
      } else { // collinear column
        pvt.shiftToEnd(j);
        pvt.setSize(pvt.size() - 1);
        k0--;
        j--;
      }
    }

    // initial QR transformation
    for (int j = k0; j < Math.min(pvt.size(), m); j++) {
      if (sum2(p[j], ks, m - 1, true) > TINY) {
        steplsqr(b, pvt, ks, j, true);
        ks++;
      } else { // collinear column
        pvt.shiftToEnd(j);
        pvt.setSize(pvt.size() - 1);
        j--;
      }
    }

    b.m = m = ks; // reset number of rows
    pvt.setSize(ks);
  }

  /**
   * QR transformation for a least squares problem <br/>
   * A x = b <br/>
   * implicitly both A and b are transformed. pvt.size() is the psuedo-rank of
   * A.
   * 
   * @param b PaceMatrix b
   * @param pvt pivoting vector
   * @param k0 the first k0 columns of A (indexed by pvt) are pre-chosen. (But
   *          subject to rank examination.)
   * 
   *          For example, the constant term may be reserved, in which case k0 =
   *          1.
   **/
  public void lsqrSelection(PaceMatrix b, IntVector pvt, int k0) {
    int numObs = m; // number of instances
    int numXs = pvt.size();

    lsqr(b, pvt, k0);

    if (numXs > 200 || numXs > numObs) { // too many columns.
      forward(b, pvt, k0);
    }
    backward(b, pvt, pvt.size(), k0);
  }

  /**
   * Sets all diagonal elements to be positive (or nonnegative) without changing
   * the least squares solution
   * 
   * @param Y the response
   * @param pvt the pivoted column index
   */
  public void positiveDiagonal(PaceMatrix Y, IntVector pvt) {

    int[] p = pvt.getArray();
    for (int i = 0; i < pvt.size(); i++) {
      if (A[i][p[i]] < 0.0) {
        for (int j = i; j < pvt.size(); j++) {
          A[i][p[j]] = -A[i][p[j]];
        }
        Y.A[i][0] = -Y.A[i][0];
      }
    }
  }

  /**
   * Stepwise least squares QR-decomposition of the problem A x = b
   * 
   * @param b PaceMatrix b
   * @param pvt pivoting vector
   * @param ks number of transformed columns
   * @param j pvt[j], the column to adjoin or delete
   * @param adjoin to adjoin if true; otherwise, to delete
   */
  public void steplsqr(PaceMatrix b, IntVector pvt, int ks, int j,
    boolean adjoin) {
    final int kp = pvt.size(); // number of columns under consideration
    int[] p = pvt.getArray();

    if (adjoin) { // adjoining
      int pj = p[j];
      pvt.swap(ks, j);
      double dq[] = h1(pj, ks);
      int pk;
      for (int k = ks + 1; k < kp; k++) {
        pk = p[k];
        h2(pj, ks, dq[1], this, pk);
      }
      h2(pj, ks, dq[1], b, 0); // for matrix. ???
      A[ks][pj] = dq[0];
      for (int k = ks + 1; k < m; k++) {
        A[k][pj] = 0;
      }
    } else { // removing
      int pj = p[j];
      for (int i = j; i < ks - 1; i++) {
        p[i] = p[i + 1];
      }
      p[ks - 1] = pj;
      double[] cs;
      for (int i = j; i < ks - 1; i++) {
        cs = g1(A[i][p[i]], A[i + 1][p[i]]);
        for (int l = i; l < kp; l++) {
          g2(cs, i, i + 1, p[l]);
        }
        for (int l = 0; l < b.n; l++) {
          b.g2(cs, i, i + 1, l);
        }
      }
    }
  }

  /**
   * Solves upper-triangular equation <br/>
   * R x = b <br/>
   * On output, the solution is stored in b
   * 
   * @param b the response
   * @param pvt the pivoting vector
   * @param kp the number of the first columns involved
   */
  public void rsolve(PaceMatrix b, IntVector pvt, int kp) {
    if (kp == 0) {
      b.m = 0;
    }
    int i, j, k;
    int[] p = pvt.getArray();
    double s;
    double[][] ba = b.getArray();
    for (k = 0; k < b.n; k++) {
      ba[kp - 1][k] /= A[kp - 1][p[kp - 1]];
      for (i = kp - 2; i >= 0; i--) {
        s = 0;
        for (j = i + 1; j < kp; j++) {
          s += A[i][p[j]] * ba[j][k];
        }
        ba[i][k] -= s;
        ba[i][k] /= A[i][p[i]];
      }
    }
    b.m = kp;
  }

  /**
   * Returns a new matrix which binds two matrices together with rows.
   * 
   * @param b the second matrix
   * @return the combined matrix
   */
  public PaceMatrix rbind(PaceMatrix b) {
    if (n != b.n) {
      throw new IllegalArgumentException("unequal numbers of rows.");
    }
    PaceMatrix c = new PaceMatrix(m + b.m, n);
    c.setMatrix(0, m - 1, 0, n - 1, this);
    c.setMatrix(m, m + b.m - 1, 0, n - 1, b);
    return c;
  }

  /**
   * Returns a new matrix which binds two matrices with columns.
   * 
   * @param b the second matrix
   * @return the combined matrix
   */
  public PaceMatrix cbind(PaceMatrix b) {
    if (m != b.m) {
      throw new IllegalArgumentException("unequal numbers of rows: " + m
        + " and " + b.m);
    }
    PaceMatrix c = new PaceMatrix(m, n + b.n);
    c.setMatrix(0, m - 1, 0, n - 1, this);
    c.setMatrix(0, m - 1, n, n + b.n - 1, b);
    return c;
  }

  /**
   * Solves the nonnegative linear squares problem. That is,
   * <p>
   * <center> min || A x - b||, subject to x >= 0. </center>
   * <p>
   * 
   * For algorithm, refer to P161, Chapter 23 of C. L. Lawson and R. J. Hanson
   * (1974). "Solving Least Squares Problems". Prentice-Hall.
   * 
   * @param b the response
   * @param pvt vector storing pivoting column indices
   * @return solution
   */
  public DoubleVector nnls(PaceMatrix b, IntVector pvt) {
    int j, t, counter = 0, jm = -1, n = pvt.size();
    double ma, max, alpha, wj;
    int[] p = pvt.getArray();
    DoubleVector x = new DoubleVector(n);
    double[] xA = x.getArray();
    PaceMatrix z = new PaceMatrix(n, 1);
    PaceMatrix bt;

    // step 1
    int kp = 0; // #variables in the positive set P
    while (true) { // step 2
      if (++counter > 3 * n) {
        throw new RuntimeException("Does not converge");
      }
      t = -1;
      max = 0.0;
      bt = new PaceMatrix(b.transpose());
      for (j = kp; j <= n - 1; j++) { // W = A' (b - A x)
        wj = bt.times(0, kp, m - 1, this, p[j]);
        if (wj > max) { // step 4
          max = wj;
          t = j;
        }
      }

      // step 3
      if (t == -1) {
        break; // optimum achieved
      }

      // step 5
      pvt.swap(kp, t); // move variable from set Z to set P
      kp++;
      xA[kp - 1] = 0;
      steplsqr(b, pvt, kp - 1, kp - 1, true);
      // step 6
      ma = 0;
      while (ma < 1.5) {
        for (j = 0; j <= kp - 1; j++) {
          z.A[j][0] = b.A[j][0];
        }
        rsolve(z, pvt, kp);
        ma = 2;
        jm = -1;
        for (j = 0; j <= kp - 1; j++) { // step 7, 8 and 9
          if (z.A[j][0] <= 0.0) { // alpha always between 0 and 1
            alpha = xA[j] / (xA[j] - z.A[j][0]);
            if (alpha < ma) {
              ma = alpha;
              jm = j;
            }
          }
        }
        if (ma > 1.5) {
          for (j = 0; j <= kp - 1; j++) {
            xA[j] = z.A[j][0]; // step 7
          }
        } else {
          for (j = kp - 1; j >= 0; j--) { // step 10
            // Modified to avoid round-off error (which seemingly
            // can cause infinite loop).
            if (j == jm) { // step 11
              xA[j] = 0.0;
              steplsqr(b, pvt, kp, j, false);
              kp--; // move variable from set P to set Z
            } else {
              xA[j] += ma * (z.A[j][0] - xA[j]);
            }
          }
        }
      }
    }
    x.setSize(kp);
    pvt.setSize(kp);
    return x;
  }

  /**
   * Solves the nonnegative least squares problem with equality constraint. That
   * is,
   * <p>
   * <center> min ||A x - b||, subject to x >= 0 and c x = d. </center>
   * <p>
   * 
   * @param b the response
   * @param c coeficients of equality constraints
   * @param d constants of equality constraints
   * @param pvt vector storing pivoting column indices
   * @return the solution
   */
  public DoubleVector nnlse(PaceMatrix b, PaceMatrix c, PaceMatrix d,
    IntVector pvt) {
    double eps = 1e-10 * Math.max(c.maxAbs(), d.maxAbs())
      / Math.max(maxAbs(), b.maxAbs());

    PaceMatrix e = c.rbind(new PaceMatrix(times(eps)));
    PaceMatrix f = d.rbind(new PaceMatrix(b.times(eps)));

    return e.nnls(f, pvt);
  }

  /**
   * Solves the nonnegative least squares problem with equality constraint. That
   * is,
   * <p>
   * <center> min ||A x - b||, subject to x >= 0 and || x || = 1. </center>
   * <p>
   * 
   * @param b the response
   * @param pvt vector storing pivoting column indices
   * @return the solution
   */
  public DoubleVector nnlse1(PaceMatrix b, IntVector pvt) {
    PaceMatrix c = new PaceMatrix(1, n, 1);
    PaceMatrix d = new PaceMatrix(1, b.n, 1);

    return nnlse(b, c, d, pvt);
  }

  /**
   * Generate matrix with standard-normally distributed random elements
   * 
   * @param m Number of rows.
   * @param n Number of colums.
   * @return An m-by-n matrix with random elements.
   */
  public static Matrix randomNormal(int m, int n) {
    Random random = new Random();

    Matrix A = new Matrix(m, n);
    double[][] X = A.getArray();
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        X[i][j] = random.nextGaussian();
      }
    }
    return A;
  }

  /**
   * Returns the revision string.
   * 
   * @return the revision
   */
  @Override
  public String getRevision() {
    return RevisionUtils.extract("$Revision$");
  }

  /**
   * for testing only
   * 
   * @param args the commandline arguments - ignored
   */
  public static void main(String args[]) {
    System.out.println("================================================"
      + "===========");
    System.out.println("To test the pace estimators of linear model\n"
      + "coefficients.\n");

    double sd = 2; // standard deviation of the random error term
    int n = 200; // total number of observations
    double beta0 = 100; // intercept
    int k1 = 20; // number of coefficients of the first cluster
    double beta1 = 0; // coefficient value of the first cluster
    int k2 = 20; // number of coefficients of the second cluster
    double beta2 = 5; // coefficient value of the second cluster
    int k = 1 + k1 + k2;

    DoubleVector beta = new DoubleVector(1 + k1 + k2);
    beta.set(0, beta0);
    beta.set(1, k1, beta1);
    beta.set(k1 + 1, k1 + k2, beta2);

    System.out.println("The data set contains " + n + " observations plus "
      + (k1 + k2) + " variables.\n\nThe coefficients of the true model"
      + " are:\n\n" + beta);

    System.out.println("\nThe standard deviation of the error term is " + sd);

    System.out.println("==============================================="
      + "============");

    PaceMatrix X = new PaceMatrix(n, k1 + k2 + 1);
    X.setMatrix(0, n - 1, 0, 0, 1);
    X.setMatrix(0, n - 1, 1, k1 + k2, random(n, k1 + k2));

    PaceMatrix Y = new PaceMatrix(X.times(new PaceMatrix(beta)).plusEquals(
      randomNormal(n, 1).times(sd)));

    IntVector pvt = IntVector.seq(0, k1 + k2);

    /*
     * System.out.println( "The OLS estimate (by jama.Matrix.solve()) is:\n\n" +
     * (new PaceMatrix(X.solve(Y))).getColumn(0) );
     */

    X.lsqrSelection(Y, pvt, 1);
    X.positiveDiagonal(Y, pvt);

    PaceMatrix sol = (PaceMatrix) Y.clone();
    X.rsolve(sol, pvt, pvt.size());
    DoubleVector betaHat = sol.getColumn(0).unpivoting(pvt, k);
    System.out
      .println("\nThe OLS estimate (through lsqr()) is: \n\n" + betaHat);

    System.out
      .println("\nQuadratic loss of the OLS estimate (||X b - X bHat||^2) = "
        + (new PaceMatrix(X.times(new PaceMatrix(beta.minus(betaHat)))))
          .getColumn(0).sum2());

    System.out.println("============================================="
      + "==============");
    System.out.println("             *** Pace estimation *** \n");
    DoubleVector r = Y.getColumn(pvt.size(), n - 1, 0);
    double sde = Math.sqrt(r.sum2() / r.size());

    System.out.println("Estimated standard deviation = " + sde);

    DoubleVector aHat = Y.getColumn(0, pvt.size() - 1, 0).times(1. / sde);
    System.out.println("\naHat = \n" + aHat);

    System.out.println("\n========= Based on chi-square mixture ============");

    ChisqMixture d2 = new ChisqMixture();
    int method = MixtureDistribution.NNMMethod;
    DoubleVector AHat = aHat.square();
    d2.fit(AHat, method);
    System.out.println("\nEstimated mixing distribution is:\n" + d2);

    DoubleVector ATilde = d2.pace2(AHat);
    DoubleVector aTilde = ATilde.sqrt().times(aHat.sign());
    PaceMatrix YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
    X.rsolve(YTilde, pvt, pvt.size());
    DoubleVector betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
    System.out.println("\nThe pace2 estimate of coefficients = \n" + betaTilde);
    System.out.println("Quadratic loss = "
      + (new PaceMatrix(X.times(new PaceMatrix(beta.minus(betaTilde)))))
        .getColumn(0).sum2());

    ATilde = d2.pace4(AHat);
    aTilde = ATilde.sqrt().times(aHat.sign());
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
    X.rsolve(YTilde, pvt, pvt.size());
    betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
    System.out.println("\nThe pace4 estimate of coefficients = \n" + betaTilde);
    System.out.println("Quadratic loss = "
      + (new PaceMatrix(X.times(new PaceMatrix(beta.minus(betaTilde)))))
        .getColumn(0).sum2());

    ATilde = d2.pace6(AHat);
    aTilde = ATilde.sqrt().times(aHat.sign());
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
    X.rsolve(YTilde, pvt, pvt.size());
    betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
    System.out.println("\nThe pace6 estimate of coefficients = \n" + betaTilde);
    System.out.println("Quadratic loss = "
      + (new PaceMatrix(X.times(new PaceMatrix(beta.minus(betaTilde)))))
        .getColumn(0).sum2());

    System.out.println("\n========= Based on normal mixture ============");

    NormalMixture d = new NormalMixture();
    d.fit(aHat, method);
    System.out.println("\nEstimated mixing distribution is:\n" + d);

    aTilde = d.nestedEstimate(aHat);
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
    X.rsolve(YTilde, pvt, pvt.size());
    betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
    System.out.println("The nested estimate of coefficients = \n" + betaTilde);
    System.out.println("Quadratic loss = "
      + (new PaceMatrix(X.times(new PaceMatrix(beta.minus(betaTilde)))))
        .getColumn(0).sum2());

    aTilde = d.subsetEstimate(aHat);
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
    X.rsolve(YTilde, pvt, pvt.size());
    betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
    System.out
      .println("\nThe subset estimate of coefficients = \n" + betaTilde);
    System.out.println("Quadratic loss = "
      + (new PaceMatrix(X.times(new PaceMatrix(beta.minus(betaTilde)))))
        .getColumn(0).sum2());

    aTilde = d.empiricalBayesEstimate(aHat);
    YTilde = new PaceMatrix((new PaceMatrix(aTilde)).times(sde));
    X.rsolve(YTilde, pvt, pvt.size());
    betaTilde = YTilde.getColumn(0).unpivoting(pvt, k);
    System.out.println("\nThe empirical Bayes estimate of coefficients = \n"
      + betaTilde);

    System.out.println("Quadratic loss = "
      + (new PaceMatrix(X.times(new PaceMatrix(beta.minus(betaTilde)))))
        .getColumn(0).sum2());

  }
}
